require(data.table)


rm(list=ls())


# This script measures the curvature and angle of the stickleback's body during the three point bending test.


dir <- "All 3P Bending Curves"
setwd(dir)


#Read the names of the .csv files in the current folder
fileNames <- list.files(pattern="*.csv", all.files=TRUE, no..=TRUE)
#Open all the files
options(datatable.fread.datatable=FALSE)
list2env(
  lapply(setNames(fileNames, make.names(gsub("*.csv$", "", fileNames))), 
         fread), envir = .GlobalEnv)


#Make a data frame to hold final data
Final3PBendCurves <- data.frame(Treatment=character(), Animal= double(), Cal=double(),
                                curve0.5=double(),curve1.0=double(),curve1.5=double(),curve2.0=double(),
                                curve2.5=double(),curve3.0=double(),curve3.5=double())

for(i in 1:length(fileNames)){
  nameNow <- fileNames[i]
  #Get rid of ".csv" in the filename
  nameNow2 <- substr(nameNow,1,nchar(nameNow)-4)
  #Open the file
  dataNow <- eval(parse(text= nameNow2))
  
  names(dataNow) <- c('X','Y')
  
  #Extract video data from filename, write it to final data frame
  splitStrings <- strsplit(nameNow2,fixed=FALSE, split = '_')
  
  treatNow <- splitStrings[[1]][2]
  Final3PBendCurves[i,1] <- treatNow
  
  aNow <- as.integer(substr(splitStrings[[1]][3],2,3))
  Final3PBendCurves[i,2] <- aNow
  
  
  #First get calibration values from each image and average them to get the overall calibration
  seq1 <- seq(1,31,5)
  seq2 <- seq(2,32,5)
  seqAll <- sort(c(seq1,seq2))
  
  calVals <- vector()
  for (cali in 1:7) {
    #All calibration distances are 1 cm. So calNow is in units pixels/cm
    calNow <- sqrt((dataNow[seqAll[2*cali-1],1]-dataNow[seqAll[2*cali],1])^2+(dataNow[seqAll[2*cali-1],2]-dataNow[seqAll[2*cali],2])^2)
    calVals[cali] <- calNow
  }
  avgCal <- mean(calVals)
  
  Final3PBendCurves[i,3] <- avgCal
  
  #adjust all X and Y values to be in cm for the curvature calculation
  dataNowCal <- dataNow/avgCal
  
  #Now we go through each bending point (0.5-3.5cm by increments of 0.5cm) and find the curvature and angle.
  
  for (k in 1:7) {
    
    x1 <- dataNowCal[5*k-2,1]
    x2 <- dataNowCal[5*k-1,1]
    x3 <- dataNowCal[5*k,1]
    
    y1 <- dataNowCal[5*k-2,2]
    y2 <- dataNowCal[5*k-1,2]
    y3 <- dataNowCal[5*k,2]
    
    
    #Curvature first
    length1 <- sqrt((x2-x1)^2+(y2-y1)^2)
    length2 <- sqrt((x3-x2)^2+(y3-y2)^2)
    length3 <- sqrt((x1-x3)^2+(y1-y3)^2)
    
    areaTri <- abs((x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2)
    
    #if area=0 then the points are co-linear and the curvature will be infinite, so we set it to NaN instead. Otherwise, we calculate curvature
    if (areaTri==0) {
      mCurveNow <- 'NaN'
    } else {
      mCurveNow <- (4*areaTri)/(length1*length2*length3)
    }
    Final3PBendCurves[i,3+k] <- mCurveNow
    
    
    
    #calculate linear distance between razors
    Final3PBendCurves[i,] razDist <- 
  }
  
}

write.csv(Final3PBendCurves, file="final3PBendCurvatureData.csv")  
  
  
  
  
  
  
  